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INTRODUCTION 


The  tschnica.!  content  of  the  resea.rch  carried  out  under  the  aegis  ot  this  contract  is 

covered  under  two  essentially  distinct  topics; 

A.  Inversion  problems  associated  with  surtace  acoustic  waves  (SAW) 

B.  Inversion  problems  associated  with  laser  doppler  velocimetry  (LDV);  signal-to-noise 

issues  occurring  in  LDV. 

These  two  topics  will  be  discussed  separately. 

We  should  point  out  that  no  attempt  was  made  to  provide  detailed  numerical  results  as 
has  been  my  custom  for  the  last  thirty  years.  The  reasons  are  two-told:  tirstlyy  these  two  topics 
do  not  lend  themselves  to  a  numerical  description  in  terms  of  a  tew  parameters  (especially  LDV) 
so  that  really  detailed  parameter  calculations  must  be  carried  (it  is  my  feeling  that  a  limited 
parameter  series  of  calculations  will  probably  lead  to  such  an  incomplete  numerical  description 
that  false  "rule-of-thumb"  conclusions  may  be  drawn);  secondly,  the  mathematical  desciiption  of 
both  topics  leads  to  nonlinear  mathematical  operations  that  are  just  too  complicated  (and  require 
enormous  memory)  to  be  run  on  workstations  of  even  VAXs,  we  need  massively  parallel 
computing  power  to  effect  the  numerical  calculations  such  as  a  CONNECTION  machine  or  the 
new  IBM  parallel  machine  on  MAUI.  Nevertheless  numerical  computations  were  carried  out  on 
a  VAX  to  illustrate  some  aspects  of  topic  A. 


NUMERICAL  INVERSION  OF  NONLINEAR  INTEGRAL  EQUATIONS  OF  TFIE  FIRST 
KIND  WITFI  APPLICATIONS  TO  SURFACE  ACOUSTIC  WAVE  PROBLEMS 


The  second  topic;  inversion  of  nonlinear  integral  equations  of  the  first  kind  with 
applications  to  SAW  (surface  acoustic  waves),  has  been  treated  in  a  quite  general  manner 
concentrating  upon  the  development  of  an  algorithm  that  is  reasonably  robust  with  respect  to 
measurement  noise.  The  technique  developed  depends  upon  a  relatively  new  idea,  that  of  "trust 
regions",  in  conjunction  with  Newton-like  methods  which  involve  exact  calculations  of  both 
Jacobian  and  Hessians  of  the  least  squares  objective  function.  The  necessary  reasons  for  using 
this  approach  (as  well  as  the  mathematical)  are  worked  out  and  discussed  in  the  main  text  of  this 
section.  Rather  than  attempting  the  inversion  problem  in  the  context  ot  SAW,  I  ha\-e  decided  to 
cast  the  specifics  in  terms  of  a  formally  similar  problem;  inversion  of  the  modulus/phase  of  a 
coherently  illuminated  object  from  its  measured  diffraction  image"  as  I  have  worked  extensively 
in  the  area  of  optical  diffraction  theory  with  the  result  that  I  have  a  good  heuristic  feeling  for  the 
plausibility  of  the  correctness  of  the  "solution".  I  am  now  confident  ol  the  usefulness  of  the 
proposed  algorithm  having  run  some  numerical  inversions  (see  main  text).  It  is  now'  a 
straightforward  matter  to  transform  the  algorithm  specifics  into  the  SAW  inversion  problem  in 
accordance  w'ith  information  to  be  provided  by  Dr.  W.  Micelli. 

As  noted  in  the  last  paragraph  of  page  20,  a  constrained  version  of  the  problem  has  been 
developed  in  order  to  handle  (if  necessary)  the  non-negativity  of  the  solution.  Further,  as  noted 
at  the  end  of  page  8,  we  really  need  to  use  a  parallel  machine  (such  as  on  Maui)  if  we  are  going 
to  speed  up  the  calculations  for  the  SAW  problem.  The  calculations  presented  here  were  carried 
out  on  a  VAX,  a  truly  slow  machine  by  todays  standards. 
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ABSTRACT 

An  algorithm  for  recovering  the  modulus  and  phase  of  a  coherently  illuminated  object 
from  its  Measured  diffraction  image  is  presented.  The  algorithm  is  based  upon  the  fact 
that  both  the  Jacobian  and  Hessian  matrices  can  be  evaluated  exactly  so  that  both  slope 
and  curvature  Information  is  available.  The  inversion  problem  is  cast  as  a  noidinear  uncon¬ 
strained  optimization  problem,  and  trust  region  techniques  are  employed  for  its  solution. 
Representative  numericals  are  presented. 
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1.  INTRODUCTION 


The  problem  of  measuring  the  diffraction  image  of  an  incoherently  illuminated 
object  and  working  back  to  determine  (or  estimate)  the  object  itsell  has  been  ot  great 
interest  in  many  scientific  and  technical  areas  for  the  past  twent}-  years.  The  subject  has 
orown  to  such  an  extent  that  even  listing  the  books  and  major  review  articles  is  a  tedious 

undertaking. 

At  the  other  extreme,  we  have  the  analogous  problem  tor  a  coherently  illuminated 
object,  a  much  more  difficult  problem  because  the  relation  betw  een  object  and  dittraction 
imase  is  nonlinear,  whereas  the  corresponding  relation  tor  incoherently  illuminated 
objects  is  linear.  In  many  respects  the  coherently  illuminated  situation  is  much  older  than 
the  incoherently  illuminated  situation  in  that  light  microscopists  have  always  encountered 
such  problems.  Their  "solutions"  have  generally  been  ot  the  old-fashioned  expert 
systems  type;  they  have  encountered  over  the  years  a  variety  ot  biological  objects  and 
developed  an  empirical  expertise  in  sorting  out  situations.  In  a  sense  their  primary 
artifact  is  the  diffraction  image,  not  eh  actual  object,  as  witness  the  various  phase  contrast 
methods  [1,2].  Unlike  the  incoherent  situation,  there  is  no  guarantee  that  the  topology  of 
the  object  bears  any  resemblance  to  its  coherent  diffraction  image.  As  an  example,  see 
Fisure  1  which  shows  the  diffraction  image  of  an  edge  viewed  through  an  optical  system 

with  an  annular  aperture  [3]. 
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In  the  present  paper,  a  solution  for  the  modulus  and  phase  ot  a  coherently 
illuminated  object  is  obtained  via  a  nonlinear  regularized  minimization  algorithm.  W'e 
have  brousht  to  bear  powerful  tools  recently  developed  in  numerical  analysis  (particularly 
trust  region  considerations)  toward  the  efficient  solution  of  the  inversion  problem.  The 
present  scheme  takes  advantage  of  the  fact  that  we  can  evaluate  both  Jacobian  and 
Hessian  matrices  exactly,  thereby  allowing  use  of  slope  and  curvature  information 
explicitly.  There  is  no  need  to  make  the  usual  small  residual  approximation  of  the 
Hessian  with  its  coiwergence  limitation  in  the  presence  of  measurement  noise.  A  second 
benefit  of  knowing  the  Jacobian  and  Hessian  e.xplicitly  is  that  we  can  make  very  efficient 
use  of  the  trust  region  tests  in  determining  the  pat  to  local  minima.  There  are  a  number  of 
alsorithms  for  the  inversion  of  modulus/phase  ot  coherently  illuminated  objecti  already 
published;  see  Stark  [4]  for  an  extremely  useful  summary.  Practically  all  these  algorithms 
will  yield  a  reasonably  accurate  estimate  of  the  support  of  the  object,  and  the  present 
algorithm  is  no  exception.  However,  many  coherently  illuminated  objects  contain 
changes  in  modular/phase  over  the  surface  so  that  an  algorithm  (such  as  the  one  discussed 
here)  which  also  can  deliver  information  on  the  distribution  of  modulus/phase  over  the 
object  is  extremely  useful  for  applications  (e.g.,  biological  light  microscopy). 
Furthermore,  the  present  method  is  reasonably  robust  with  respect  to  noise  because  we 
operate  in  the  large  residual  regime. 


2.  DIFFRACTION  MODEL 


The  diffraction  image  of  a  coherently  illuminated  object,  J(x,y)  in  the  image  receiving 
plane  with  coordinates  x,  y  is  measured  over  a  square  lattice  of  points  x^,  Vn- 

^  0,±1,±2,...,±T/  (2.1) 

Vn  =  Pn 

where  ,5  is  a  numerical  constant.  It  is  assumed  that  M  is  large  enough  for 


In  what  follows,  it  will  be  convenient  to  write  (x„,,ym)  =  ^mn  as  a  column  vector  I  of 
length  m  =  (2.\/  +  1)'  using  standard  Fortran  lexiographic  ordering. 

We  assume  that  the  measured  diffraction  image  can  be  modelled  via  scalar  optical 
diffraction  theory  so  that  the  model  diffraction  image,  I(x,y).  is  given  by  the  convolution 
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object 


(2.3) 


assurning  the  isoplcinatic  condition  to  hold. 

The  coherently  illuminated  object  is  characterized  by  a  complex-valued  function 
o(x,y)  with  modulus  lo(x,y)l  and  phase  arctan[o.-(x, y)/o„(x, y)].  The  function  a(x,y) 
is  the  coherent  point-spread  function  of  the  optical  system  performing  the  imaging;  to 
within  multiplicative  factors  it  is 


a(,.v)=  Jj 

exit  pupil 


(2.4) 


S 


with  k  cLS  the  mean  wavenumber  of  the  coherent  light  and  /  the  focal  length.  The  pupil 
function  is  given  by 


(2.5) 


with  W{C,ri)  wavefront  aberration  function  over  the  exit  pupil  and  B{C,r])  the  ampli¬ 
tude  distribution  over  the  e.xit  pupil;  both  are  assumed  known  in  the  present  scenario. 


For  the  very  important  case  of  a  circular  aperture  e.xit  pupil  of  radius  a  for  which 
B{C,rj)  =  1  and  T'F((;,  r?)  =  0,  we  have 


a{x,y) 


2.7i 


¥(x- 


2  ^  1  /2 
y  )  ' 


(2.6) 


Since  in  most  applications,  the  exit  pupil  is  circular;  under  these  circumstances  it  is 
convenient  to  employ  normalized  coordinates; 


kci 

u  =  -j-x, 

c 


p  = 


a 


ka 

V  =  -yV 
V 


(2.7: 


Consequently  Eqs.  (2.3),  (2.4),  and  (2.6)  become 


I{u,v)  = 


J  j  a{u  -  u',v  —  v')o(u',v')  du'dv' 


object 


a 


(2.3a) 


a{u,v) 


j  J  ^(p,5)e’(“P+''’^  dpdq 


a(u, v) 


2Ji((u^  +  t'")^'^^) 

(u2+ ^2)1/2 


(2.4a) 


(2.6a) 


Leaving  aside  the  mathematical  details  for  the  moment,  let  us  consider  the  task  before 
us.  We  are  given  measured  values  of  the  diffracted  intensity  and  are  required  to  determine 
the  modulus  and  phase  of  the  object  assuming  that  the  measured  diffracted  intensity 
can  be  modeled  as  the  convolution.  Eq.  (2.3),  and  that  we  know  the  parameters  of  the 
optical  system  characterizing  the  coherent  point  spread  function,  Eq.  (2.4).  Said  object 
is  contained  in  the  integrand  of  a  double  integral  which  is  itself  squared.  When  viewed 
from  this  perspective,  we  should  not  be  unduly  optimistic  about  achieving  an  accurate 
solution  because  of  the  inherent  nonlinearity  and  strong  smoothing  action  of  the  double 
integral.  Irrespective  of  the  actual  inversion  method,  the  strong  smoothing  action  of  the 
integration  means  that  much  high  frequency  data  is  lost  and  cannot  be  retrieved,  even  in 
principle.  Only  a  low  frequency  version  of  the  object  can  be  obtained.  Perhaps  the  best 
way  to  consider  the  problem  is  to  interpret  it  thusly:  we  are  given  the  “answer’  (=  effect) 
in  the  form  of  a  noisy  diffraction  image  of  the  object  and  are  attempting  to  determine 
the  “question”  (=  cause),  the  coherently  illuminated  object.  The  mathematical  relation 
between  answer  and  question  is  highly  nonlinear  by  virtue  of  Eq.  (2.2),  bearing  in  mind 
that  a{u,v)  is  additionally  an  oscillating,  complex- valued  function. 


In  a  sense  this  problem  can  be  considered  as  a  two-dimensional  phase  retrieval  problem, 
see  for  various  details.  We  will  not  discuss  the  general  issues  of  ill-posed  inverse 

problems  and  refer  the  reader  to  [53^]  for  some  of  the  theoretical  aspects. 


Our  diffraction  model  image  I{u,  v)  is  to  be  determined  by  the  values  of  the  .image 
at  the  lattice  points  inside  a  square  in  the  {u',v')  plane  that  is  large  enough  to 

contain  the  image.  The  number  of  object  values,  o(tq,,u^)  —  is  taken  to  be  N. 
important  to  remember  that  both  modulus  and  phase  must  be  determined  at  each  of  these 
lattice  points.  Let  h  be  the  number  of  unknown  modulus  and  phase  values  (obviously  n 
must  be  even);  the  unknown  Ok^  are  to  be  written  as  a  column  vector  o. 


Although  we  have  run  several  dozen  inversions,  we  do  not  possess  a  sufficiently  large 
base  to  estimate  benchmark  performance.  However,  many  of  our  inversions,  regardless 
of  the  initial  guess,  took  150-1S0  iterates  to  converge.  In  a  few  cases,  convergence  was 
achieved  in  half  this  number;  nevertheless  we  also  encountered  some  cases  where  conver¬ 
gence  required  200  iterations. 

Generally  the  algorithm  performance  decomposed  into  two  stages:  (a)  the  initial  stage 
wherein  the  objective  functions  went  through  wild  gyrations  in  magnitude  as  the  trust  re¬ 
gion  subalgorithm  had  to  compute  many  new  gradient  and  Hessian  matrices  while  searching 
for  an  appropriate  direction  of  descent;  (b)  the  second  stage,  when  the  trust  region  subal- 
o-orithm  has  locked  into  a  subset  of  '‘optimal”  directions  of  descent,  the  objective  function 

O 

then  undergoes  a  monotone  decrease  to  zero  as  the  number  of  iterations  is  increased.  We 
again  warn  the  reader  that  although  convergence  is  attained,  there  is  no  guarantee  that 
the  global  minimum  has  been  achieved. 

As  a  final  remark,  we  note  that  as  the  analysis  is  heavily  dependent  upon  linear  algebra 
(matrices,  vectors);  it  would  be  of  some  interest  to  consider  doing  calculations  via  parallel 
processing  algorithms  to  speed  up  the  computations  e\en  more. 


In  this  notation  Eq.  (2.3a)  becomes  ^(o)  where  J  is  a  vector  of  length  m.  The 
inversion  problem  relates  the  unknown  (complex-valued)  object  values  o  of  the  assumed 
diffraction  model  to  the  measured  diffraction  image  I  via  the  nonlinear  relation 


I  =  J(o) 


(2.7) 


This  equation  is  to  be  interpreted  as  a  system  of  m  nonlinear  equations  in  h  unknowns. 
Enough  data  must  be  given  to  allow  some  smoothing  of  tire  measured  diffraction  image 
data  as  regards  the  diffraction  model;  consequently  we  let  m  >  n  so  that  the  nonlinear 
system  is  overdetermined.  The  problem  now  reduces  to  the  solution  oi  Eq.  (2.1)  in  some 
normed  sense. 

It  is  of  some  interest  to  contrast  the  differences  between  the  incoherently  illuminated 
and  coherently  illuminated  object  situations.  The  relation  between  object  and  immge  for 
incoherently  illuminated  objects,  again  assuming  that  the  isoplanatic  conditions  holds,  is 

U{x,y)  =  j  J  t{x-x',y-y')o,{x',y')dx'dy'  (2.S) 

object 


where 


o,{x,y)  =  intensity  distribution  over  the 
incoherently  radiating  object 


t{x,y)  <x  la(x,y)l“  =  incoherent  point-spread  function 

J;(x,  y)  =  intensity  distribution  over  the  dirtraction 
image  of  Oi{x,  y) 


Two  points  to  note:  (1)  all  functions  in  Eq.  (2.8)  are  real  and  nonnegative,  (2)  the  relation 
between  answer,  Ii{x,  y),  and  question,  Oi(x,  y),  is  linear.  On  the  other  hand,  for  the  coher¬ 
ently  illuminated  object  only  the  diffracted  intensity  is  real  and  nonnegative.  Furthermore 
the  relation  between  I{x,y)  and  o(x,y)  is  nonlinear.  Consequently  the  incoherently  illu¬ 
minated  object  scenario  is  much  easier  to  invert  than  the  coherently  illuminated  object 

scenario. 
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3.  PRELIMINARIES 


In  order  to  proceed,  we  next  discretize  the  double  integral  on  the  right-hand  side  of 


Eq.  (2.3a) 


jV  N 


Xlm 


fc,  t'n  -  Vl)o{vk,Vi) 


k  =  l  t=l 


Here  u,,  tv  are  the  quadrature  points  and  are  the  corresponding  weight  factors. 

In  the  numerical  algorithm  we  employ  for  the  inversion  (see  next  section),  we  require 
the  first  and  second  derivatives  of  with  respect  to  ot,  in  order  to  form  .Jacobian 

and  Hessian  matrices.  We  need  not  vHte  out  the  explicit  formulas  because  we  employed 
a  symbol  manipulation  program^to^vlluate  them  directly  in  the  computer  from  which 
nuiTierical  values  are  ootained  inteinalH. 
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4.  OUTLINE  OF  SOLUTION  APPROACH 

We  now  outline  clie  genevo,!  features  of  our  numerical  approach  to  the  nonlinear  phase 
retrieval  problem  via  regularized  unconstrained  minimization.  Our  version  of  the  uncon¬ 
strained  minimization  problem  is  that  of  finding  the  least  value  of  an  objective  function, 
e{o).  The  term  unconstrained  indicates  that  the  variables  o  are  not  limited  in  any  way.  In 
our  application,  we  wish  to  determine  a  global  minimum  of  e,  i.e.,  a  point  o*  satisfying 

e(o)  >  e(o*),  Vo  .  (4.1) 

Unfortunately  we  uiust  be  content  with  a  local  minimum: 

e(o)>e(o‘)  Vo,  in  a  neighborhood  of  o*  .  (4.2) 

Solution  of  the  globa.1  minimum  problem  is  far  harder  than  the  solution  of  the  locad  mini¬ 
mum  problem  which  itself  is  extremely  complicated  [16,17]. 

The  derivation  of  nearly  all  methods  for  unconstrained  minimization  is  founded  on  the 
assumption  that  the  objective  function  e(o)  can  be  appro.ximated  by  a  quadratic  function 
in  the  neighborhood  of  a  minimum  point.  Thus  methods  are  sought  which  efficiently 
minimize  quadratic  functions  in  the  hope  that  they  will  also  be  effective  on  moie  general 
functions,  at  least  in  the  neighborhood  of  a  minimum.  11/ hen  first  and  second  derivatives  of 
e(o)  are  known,  such  as  in  our  case,  then  we  can  make  use  of  both  gradient  and  curvature 
information  to  effect  a  solution.  Furthermore,  both  gradient  and  curvature,  as  governed 
by  G  and  H,  see  Eqs.  (4.4)  and  (4.5),  can  be  calculated  exactly  in  the  context  of  the 
discretized  version  of  our  problem.  Thus,  we  can  avoid  many  of  the  difficulties  cissociated 
with  situations  where  G  and  H  are  known  only  approximately.  The  Taylor  seiies  expansion 
of  the  objective  function  can  be  used  to  approximate  its  minimum  value  from  points  o  near 
to  the  minimum  o*  by  moving  in  a  direction  Ao 

e(o  -r  Ao)  %  €(o)  -h  G^'(o)(Ao)  +  ^(Ao)^H(o)(Ao) 


(4.3) 


G  and  H  are  the  gradient  vector 


and  Hessian  matrix  of  the  objective  function,  respectively 


H  = 


de 

de 

de 

aoi  ’ 

do2  ’ 

dOn 

d-e 

aoi 

doidoh 

d-e 

(4.4) 


(4.5) 


Note  that  H  is  symmetric.  Both  G  and  H  are  exact  within  the  context  of  the  discretization. 
The  strategy  is  to  determine  the  vector  Ao  of  the  movements  required  to  approximate  the 
minimum  from  the  current  point  o. 

Before  proceeding  further,  it  is  important  to  state  the  conditions  for  a  given  point  o  to 
be  an  unconstrained  strong  minimum  o’  (i.e.,  a  point  o*  for  which  the  objective  function 
increases  locally  in  all  directions.  The  first  order  necessary  condition  is  [16,17] 


G(o*)  =  0  .  (^-6) 

However,  this  condition  is  not  sufficient  as  other  types  of  minima  (stationary  points)  also 
satisfy  this  condition.  The  second  order  condition  follows  from  Eq.  (4.2)  and  is  [11,12] 

(Ao)^H(o’)(Ao)  >  0  (4-7) 

which  is  sufficient  to  ensure  that  e(o*  +  Ao)  >  e(o’).  If  e(o)  f  0,  this  implies  that  H(o’) 
is  a  positive  definite  matrix.  Therefore  the  second-order  sufficient  condition  for  a  strong 
minimum  is  that  i7(o*)  should  be  positive  definite. 

Returning  to  Eq.  (4.3),  hereafter  termed  the  local  model  of  the  object  function,  we 
initiate  a  search  for  o’  by  moving  Ao.  This  involves  an  iterative  procedure  in  that  we 


start  from  some  point  o  and  choose  in  some  fashion  a  direction  Ao  in  which  the  minimum 
is  assumed  to  lie.  This  is  repeated  until  che  minimum  is  achieved  (if  possible).  Because  we 
are  employing  a  quadratic  version  of  e(o)  then  the  local  model  of  e(o)  always  allows  us  to 
find  a  solution.  However,  the  local  model  of  e(o)  is  certainly  not  a  useful  approximation 
to  e(o)  itself  except  near  a  minimum.  It  is  an  act  of  faith  that  the  local  model  and  the 
actual  model  are  “close"  in  the  vicinity  of  the  minimum.  Simple  calculations  show 

that  at  a  minimum 

Ao  = -[H(o)]-'G(o)  .  (4.8) 

This  eciuation  represents  the  appropriate  direction  Ao  to  take  to  the  minimum  o  ,  based 
upon  information  at  o.  This  equation  is  fundamental  to  all  second-order  minimization 
algorithms  (i.e.,  algorithms  employing  both  G  and  H).  However,  on  the  actual  objective 
function  surface,  such  as  we  are  faced  wdth,  the  local  model  of  e(o)  is  only  accurate  in 
the  immediate  vicinity  of  the  minimum..  This  means  that  H  may  not  be  positive  definite, 
as  required  by  Eq.  (4.7).  since  it  is  evaluated  at  points  other  than  the  minimum.  This 
situation  is  most  likely  to  occur  at  some  distance  from  the  minimum  (i.e.,  for  initial  iter¬ 
ations)  since  at  a  point  close  to  the  minimum  all  sufficiently  differentiable  functions  tend 
to  behave  as  a  quadratic  function  as  their  third  and  higher  order  derivatives  in  the  Taylor 
series  become  negligible.  We  will  return  to  the  necessity  of  keeping  H  positive  definite 

verv  shortly. 

The  classic  approach  to  limiting  step  size  during  the  iteration  is  a  line  search 
In  line  search  tactics  we  compute  a  descent  direction  and  subject  this  direction  (which  is 
generally  not  towmrd  the  unconstrained  minimum)  to  a  minimization  procedure  of  which 
details  can  be  found  in  the  above  references.  Should  the  descent  direction  satisfy  these 
criteria,  this  iteration  is  then  terminated.  We  then  repeat  the  process,  etc.  The  difficulty 
of  practical  implementation  is  two-fold.  First,  such  calculations  are  prohibitively  expen¬ 
sive  and  the  minimization  criteria  are  therefore  only  approximate  to  sa.ve  computer  time, 


vet  thev  must  be  made  precise  enough  to  ensure  reasonably  quick  convergence.  Second, 
fulfilling  these  contradictory  goals  is  something  of  a  black  art  and  the  programming  effort 
is  very  substantial,  often  occupying  up  to  two-thirds  of  the  coding  for  the  entire  optimiza¬ 
tion.  For  small  problems,  such  as  the  slit  aperture,  line  search  methods  are  very  useful 
and  were  employed  along  with  regularized  singular  value  decomposition  in  [IS]. 

Although  it  is  possible  to  use  line  search  methodology  for  the  2-D  aperture,  we  have 

chosen  to  employ  a  relatively  new  method  which  offers  many  advantages  over  the  line  search 

j  M  1  Co  ^ 

_  t.nrtir  A  particular  advantage  of  the 


methodology,  the  trust  region  tactic  A  particular  advantage  ot  the 

trust  region  approach  is  in  the  very  strong  global  convergence  properties  which  hold  with 
no  significant  restrictions  on  the  class  of  problems  to  which  they  apply  making  it  especially 
useful  for  the  nonlinear  minimization  of  phase  retrieval  in  two  dimensions.  In  addition,  the 
trust  region  philosophy  requires  less  computation  than  does  the  line  search  philosophy  in 
terms  of  gradient  and  Hessian,  but  more  computations  have  to  be  performed  on  the  local 

itiocIgI  of  th-C  objGcti\  c  ^unction. 

Suppose  we  are  at  point  o^-'),  after  the  ^--th  iteration.  Now  let  us  assume  that  there 
is  a  region  which  we  take  to  be  in  the  shape  of  a  sphere  of  radius 

^  {o:  llo-o^'-^ll  <  (4.9) 

in  which  the  local  model  of  e(o),  given  by  the  Taylor  series  Eq.  (4.2)  agrees  with  the  actual 
objective  function  in  some  sense.  One’s  obvious  choice  is  to  let 


(4.10) 


where  the  correction  minimizes  the  local  model  (A)  for  all  values  of  o^  )  -f  A  in  A^ 
An  iteration  step  in  o  is  to  be  restricted  by  the  region  of  validity  of  the  Taylor  series.  Thus 
we  compute  the  gradient  and  Hessian,  G  and  li,  appropriate  to  o^*-^  and  minimize  the 


local  model  of  €(o)  in  order  to  determine  the  radius  of  the  trust  region.  In  formal 
language,  we  seek  the  solution  of  the  pioblem. 

minimize  (local  model)  subject  to  constraint  IIAII2  <  .  (4.11) 


Before  we  can  solve  this  subproblem  it  is  necessary  to  choose  some  reasonable  criterion 
on  the  radius  of  the  trust  region.  Obviously,  the  criterion  should  not  present  undue 
restrictions,  so  that  should  be  as  large  as  possible  subject  to  some  agreed-upon  idea 

as  to  the  degree  of  closeness  of  the  local  model  of  e(o)  and  e(o).  One  way  is  to  define  the 


quality  coefncient 


rk  = 


e(oh^4  ^  _  e(o’^^'^)_ 

e,(oO-)  +  A(^-))-e,j(oO-)) 


(4.12) 


where  the  denominator  is  the  predicted  reduction  of  the  local  model  of  e,  now  call  it  e,, 
while  the  numerator  represents  the  reduction  in  the  actual  objective  function.  The  closer 
rjt  is  to  unity,  the  better  the  agreement  between  e  and  e^.  As  a  stopping  criterion  stop  if 
rk  >  -8  and  go  to  the  next  iteration.  If  rjt  <  .8  reduce  h‘'^^  and  repeat  the  calculation. 
The  literature  contains  other  stopping  criteria  but  the  one  quoted  above  seems  adequate 
for  our  purposes.  See  the  quoted  references  for  more  details. 


Thus  far  we  have  outlined  the  general  features  of  the  inversion  problem,  estimating 
o  via  minimization  of  an  objective  function  e,  as  yet  unspecified.  We  feel  that  there  are 
two  objective  functions  of  possible  interest.  The  first  is  the  usual  least  squares  {L2  norm) 

objective  function 

m 

=  ,  (4.13) 

In  this  version  of  the  problem,  we  consider  the  estimation  of  o  via  an  (unconstrained) 
nonlinear  least  squares  minimization.  A  second  objective  function  is  for  an  Lx  norm 

m 

e(o)  =  y]  \Ii  -  ^i)l 
i=:l 


1  (d 


(4.14) 


Consequently  we  have  allowed  considerable  flexibility  in  the  e  minimization  algorithm  so 
as  to  accommodate  both  objective  functions  as  well  as  others  that  may  arise.  For  the 
purposes  of  the  present  paper,  we  confine  the  discussion  to  the  L2  norm  aspects. 


Upon  defining  the  vector 


#(o)=  X~I 


(4.15) 


,vhich  is  of  length  m.  we  can  rewrite  Eq.  (4.13)  as 


e(o)  ==  $^(o)$(o)  . 


(4.16) 


The  elements  of  G  and  H  can  be  expressed  directly  in  terms  of  the  elements  #  by 
introducing  an  ancillary  matrix  J  given  by 


(4.17) 


do- 

■  m 

•  ^ 


This  matrix  is  generally  rectangular.  Thus 


i=i 


(4.18) 


G  =  2J^  $  . 


(4.19) 


The  corresponding  elements  of  H  are 


d~e  dSi  d(f>i  ^ 

dokdoi  ""  “  dok  doi  dokdo( 


(4.20) 


n 


with  k,  £  =  1,2, . . . ,  h.  The  first  term  on  the  right-hand  side  is  2J'^J  and  we  will  call  the 
second  term  S;  thus 

H  =  2J^J-hS  (4-21) 

There  is  no  guarantee  that  H  will  always  be  positive  definite  through  the  calculations; 
in  fact  H  will  become  almost  singular  as  we  approach  the  strong  minimum.  To  avoid 
this  state  of  affairs,  we  consider  a  regularized  version  of  H.  In  the  regularized  version, 
we  replace  H  by  H  -f  al  where  q,  the  regularizing  parameter,  is  a  small  nonnegative  real 
number.  This  procedure  can  remove  H  from  near  singularity  and  for  sufficient  large  q 
restore  the  positive  definite  character  of  the  Hessian.  In  linear  problem.s,  there  is  a  well 
established  method  to  estimate  q  [m].  In  the  nonlinear  problem  we  employed  q  in  the 
rano-e  0.01  <  ?  <  0.05.  There  was  not  much  difference  in  the  final  answers  as  long  as 
g  >  0;  but  setting  g  =  0  caused  considerable  numerical  instability  as  expected. 


5.  A  NUMERICAL  EXAMPLE 


It  is  not  our  intention  to  present  a  catalog  of  numerical  results  at  this  time,  yet  we  wish 
to  discuss  the  numerical  workings  of  the  algorithm  in  the  presence  of  "measurement”  noise. 
To  this  end  we  considered  a  known  object  o{u,  v)  and  from  it  calculated  the  diffraction 
image  I{u,v)  from  Eq.  (2.3a)  using  the  point-spread  function  a{u,v)  corresponding  to  an 
in-focus,  aberration-free  optical  system  as  given  by  Eq.  (2.6a).  The  measured  diffraction 
image  I(u,  u)  was  taken  to  be  given  by 

I[u,  v)  =  [1  -r  Sij{u,  v)]T{u,v)  (b-l) 

where  6  is  a  positive  constant  less  than  unity  and  i.iium,  Vn)  is  a  random  variable  uniformly 
distributed  over  (-1,  -i-l): 

m  \  <  1 

=  0  ,  \l-i\  >  1  (5-2) 

Values  of  6  used  in  the  present  calculations,  such  as  6  =  0.04  are  described  loosely  as 
4%  a  noise.  Note  that  the  noise  we  are  introducing  is  intensity  dependent  noise;  it  is  not 
intended  to  faithfully  simulate  actual  detector  noise  but  rather  to  mimic  such  noise  for  the 
purpose  of  testing  the  robustness  of  the  inversion  algorithm. 

We  have  chosen  the  following  lone  object  to  illustrate  the  calculations. 

o(ti,i;)  =  0  ,  ,  —  CO  <  V  <  —15 

= -[3-berf(t;)]  ,  -  15  <  v  <  15  (5.3) 

8 

=:0,  —  15<U<CO 

where  erf(u)  is  the  error  function.  There  are  two  reasons  for  choosing  this  object.  The 
first  reason  is  that  the  modulus  of  the  object  exhibits  a  spatial  variation.  One  of  the 


(as  yet  unstated)  requirements  on  the  algorithm  is  that  it  oe  able  to  reco'ver  reasonably 
accurate  estimates  of  the  object  photometry  in  addition  to  estimating  its  size;  as  noted  in 
the  introduction  almost  any  of  the  currently  available  algorithms  can  perform  this;  very 
few  seem  to  be  able  to  provide  reasonably  accurate  estimates  of  photometry.  The  second 
reason  for  choosing  this  object  is  that  it  is  very  small  in  width  being  slightly  less  than  Airy 
disk  radius  across  (recall  that  an  Airy  disk  radius  «  3.S2).  These  two  constraints  present 
a  severe  test  of  the  algorithm. 

Calculations  were  run  for  the  completely  unrealistic  case  of  the  noise- free,  measured 
diffraction  image  of  the  object  given  by  Eq.  (5.3).  The  inversion  algorithm  was  able  to 
return  reasonable  answers  as  compared  to  the  true  values,  ^\e  will  not  quote  any  of  these 
results  and  go  to  the  “noisy”  measurement  situation. 

In  Figure  2  we  show  a  sample  realization  of  the  reconstruction  in  the  presence  of  2% 
measurement  noise.  Xote  the  presence  of  negative  values  oi.  the  modulus  at  the  edges  of 
the  object;  however  the  photometry  of  the  reconstruction  follows  the  true  object  'veiv  well 
except  at  the  edges.  Two  sample  realizations  of  the  reconstruction  of  the  modulus  of  the 
object,  in  the  presence  of  measurement  noise  are  shown  in  Figs.  3  and  4.  Again  note 
the  unphvsical  negative  values  of  the  modulus  at  the  edges  of  the  object.  Considering 
the  fact  that  the  object  is  so  small,  the  modulus  is  in  reasonable  agreement  with  the  true 
object  in  spite  of  measurement  noise. 

The  troublesome  feature  of  the  present  inversion  algorithm  is  obviously  the  occurrence 
of  negative  modulus  values  at  the  edges  of  the  object.  The  edges  are  somewhat  unrealistic 
because  they  are  of  infinite  slope  and  the  algoritnm  tends  to  overshoot  in  the  ma.nner  of  a 
Gibbs  phenomenon.  \\'e  could,  of  course,  consider  objects  with  finite  slopes  at  the  edges 
to  minimize  the  overshoots  and  undershoots;  instead  we  have  decided  to  develop  a  vari¬ 
ant  of  the  Inversion  algorithm  using  constrained  minimization  to  surmount  this  annoying 
problem.  Details  will  be  reported  in  the  near  future. 


^0 


6.  SUMMARY 


Given  the  rather  discursive  presentation,  we  feel  it  is  useful  to  summarize  the  essentials 
of  our  approach  to  the  problem.  We  consider  the  problem  of  determining  the  modulus  and 
phase  of  the  object  at  the  lattice  points  as  one  of  unconstrained  minimization.  A  local 
(i.e.,  quadratic)  model  approach  is  used.  We  are  in  the  fortunate  position  of  employing 
both  the  Jacobian  and  Hessian,  which,  in  the  context  of  the  discretized  diffraction  model, 
can  be  evaluated  analytically!  Thus  we  make  use  of  both  slope  and  cur^-ature  information, 
other  methods  only  use  slope  information.  The  relatively  new,  and  powerful  trust  region 
algorithm  is  then  utilized  in  place  of  the  usual  line  search  algorithm.  The  trust  region 
ah-orithm  makes  full  use  of  the  local  model  and  takes  into  account  local  validity;  moreover 
the  algorithm  exhibirs  global  convergence  properties. 
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FIGURE  LEGENDS 


Fig.  1  Distribution  of  illuminence  in  the  diffraction  image  of  a  coherently  illuminated  opaque 
edge:  o(c)  =  0,  u  <  0  and  o{v)  =  1  for  v  >  0  viewed  through  an  annular  aperture 
of  obscuration  radius  £  =  0.05.  The  solid  line  is  for  the  in-focus  situation,  the  dotted 
line  is  for  a  half-wave  of  defocus.  Taken  from  Reference  #3. 

Fig.  2  Sample  realization  of  reconstruction  of  modulus  and  phase  of  coherently  illuminated 
object  (dashed  lines)  in  the  presence  of  2%  measurement  noise  . 

Fio-  3  Sample  realization  of  reconstruction  of  modulus  of  coherently  illuminated  object 
(dashed  line)  in  the  presence  of  ^  measurement  ’  noise  . 

Fig.  4  Sample  realization  of  reconstruction  of  modulus  of  coherently  illuminated  object 

°  2>7c)  „  . 

(dashed  line)  in  the  presence  of  4^  measurement  "noise  . 


LASER  DOPPLER  VELOCIMETRY  ISSUES 


The  subject  of  laser  doppler  velocimetry  (LDV)  has  been  of  interest  to  the  Navy  for  many 
years  and  for  a  variety  of  reasons  we  need  not  chronicle.  I  have  been  tasked  by  Dr.  William 
Statnick  to  investigate  theoretical  issues  relating  to  this  topic;  this  section  contains  some  of  my 
preliminary  investigations  into  the  subject,  of  which  there  are  three  subtopics  reported  at  this 

time. 

The  first  subtopic  is  devoted  to  the  subject  of  increasing  the  signal-to-noise  ratio  in 
multichannel  data  for  use  in  LDV.  Phis  I  hav'e  accomplished  b\  employing  a  detciministic 
version  of  the  KL-expansion,  the  details  are  sketched  out  in  the  report.  .A  second  subtopic  which 
arose  from  the  first  subtopic  is  that  of  determining  a  harmonic  signal  bulled  in  Gaussian  noise. 
While  working  on  the  deterministic  KL-expansion,  I  figured  a  way  to  extend  the  stochastic  KL- 
expansion  to  second  order  statistics.  The  work  has  been  written  for  possible  publication,  and  the 
details  are  given  here  in  the  form  of  a  paper.  The  third  subtopic  is  that  of  LDV  itself 


SIGNAL-TO-NOISE  ENHANCEMENT  IN  MULTICHANNEL 
DATA  FOR  USE  IN  LVD 

Dr.  William  Statnick  has  raised  the  issue  of  understanding  how  to  increase  signal  to  noise 
ratios  (SNR)  in  multichannel  data  for  use  in  LDV,  a  very-  important  topic  with  practical  bearing. 

I  have  been  able  to  develop  an  approach  to  the  SNR  enhancement  problem  using  the  Karhunen- 
Loeve  transform.  This  approach  will  now  be  sketched;  as  noted  in  the  introduction  numerical 
computations  were  not  attempted  because  of  the  necessity  of  having  to  use  a  parallel  machine 
such  as  the  CONTIECTION  machine  or  the  PS-2  IBM  machine  (at  MAUI).  Hopefully  if  a 
follow-on  is  granted,  then  I  intend  to  apply  for  time  on  the  MAUI  machine  to  perform  the 
numerical  calculations.  The  reason  is  quite  simple,  we  need  to  look  at  as  many  of  200  channels 
in  order  to  be  realistic,  serial  machines  are  simply  too  slow.  I  doubt  that  a  small  number  of 
channels  (say  10  or  20)  is  of  any  real  importance. 

The  basic  idea  behind  my  use  of  the  KL  expansion  (or  transform)  is  that  it  is  the  optimal 
way  to  extract  coherent  information  from  multichannel  input  data  (such  as  is  needed  in  some 
LDV  scenarios)  in  a  least  squares  norm  sense.  There  are  many  ways  to  derive  the  KL  expansion, 
but  the  present  approach  is  probably  the  simplest  for  our  purposes.  A  more  mathematical 
approach  using  gap  integral  equations  is  discussed  in  the  next  subsection. 

Suppose  we  have  a  set  of  n  real  signals  (we  could  consider  n  complex  signals,  but  the 
analvsis  becomes  verv  elaborate  without  any  additional  insight).  Call  these  signals  Xj(t)  where) 
=  l,2,...,n.  Define  a  transform  Xj(t)  with  an  associated  transformation  matrix  A  (as  yet 


undefined)  such  that 


XU)  =  £aijX,(t) 
^  1=1 


where  &rc  the  elements  of  mntrix  A.  We  wnnt  to  choose  the  X^(t)  such  thut  they  form  un 
orthogonal  basis,  that  is 


ti 

x,(t)  =  Zb|jX/t) 

J=1 


If  the  number  of  channels  is  very  large,  then  we  approximate  the  above  by 


m 

y,(t)  =  Zb„x,(t) 

where  m  <  n.  How  large  m  must  be  to  approximate  n  can  only  be  determined  by  numerical 
computation.  The  yx(t)  are  the  reconstructed  signals  and  by  are  the  matrix  elements  of  B,  the 
inverse  expansion  matrix. 

Our  object  is  to  reconstruct  the  Xj(t)  using  the  smallest  number  of  basis  signals  m,  ot 

A- 

course  to  some  specified  eiror.  We  choose  (given  a  fixed  m)  to  require  that  the  matiices  A  and  B 
are  such  that  they  minimize  the  least-squares  error 


yi(t)]  dt 


where  T  is  the  observation  time. 

It  can  be  shown  that  the  rows  of  the  matrix  %  consist  of  the  normalized  eigenvectors  of 
the  matrix  C 


C  is  something  like  a  covariance  matrix,  note  that  it  is  real,  symmetric,  and  positive  semi- 
defmite;  thus  it  possesses  the  decomposition 

^  A  A  A 

c  =  USU" 

A 

The  matrix  U  contains  the  normalized  eigenvectors  U;  where 

A  A  A 

CU,=XjU, 

and  S  is  a  diagonal  matrix  of  eigenvalues  (A-i,  X2,  •••,  !!•„)•  Note  that  these  eigenvalues  are 
nonnegative.  When  A=U,  the  rotated  signals  Xj(t)  form  an  n-dimensional  subspace  and  we  can 
show  that  the  least-squares  error  is  now  given  by 


e(m)  = 

j-m  +  ] 


Since  the  eigenvalues  are  arranged  in  descending  order,  it  follows  that  the  lower  order  basis 
functions  can  be  used  to  reconstruct  most  of  the  signal  (i.e.,  the  lower  order  basis  functions  are 
somewhat  like  lower  order  Fourier  basis  functions). 

We  note  the  following,  for  our  SNR  application; 

1 .  The  KL-expansion  produces  a  set  of  uncorrelated  (or  somewhat  more  precisely  in  the 
present  context,  orthogonal)  basis  functions  from  the  data  set. 

2.  The  magnitude  of  the  j-th  eigenvalue  is  a  measure  of  the  amount  of  coherent  energy 
present  in  the  j-th  basis  function  (for  the  proof  see  the  appendix). 

3.  The  fact  in  item  2  implies  that  the  reconstructed  signal  using  only  the  lower  order 
basis  functions  amounts  to  reconstruction  ot  the  coherent  energy  present  in  the  input. 


I  want  to  emphasize  that  this  version  of  the  KL-transfomi  is  deterministic  in  that  the  input 
signals  are  deterministic,  whereas  the  analysis  in  the  next  subsection  deals  with  the  stochastic 
aspects  of  the  KL-transform.  In  view  of  this  difference,  I  suggest  that  a  detailed  numerical 
solution  of  both  cases  be  carried  out  in  order  to  determine  which  case  really  applies  to  the  LDV 
scenario. 


Appendix  A 

The  energy  of  the  basis  functions  can  be  studied  thusly:  Given  the  data  vector 
x(t)  =  {Xi(t),  i  =  l,2,...,n},  compute  the  C  matrix 

A  A  ^ 

C  =  xx'' 

A 

Since  C  is  symmetric,  then 

C  =  UAU 

where  U  is  the  matrix  of  column  eigenvectors  Uj  and  A  is  the  diagonal  matrix  of  eigenvalues. 
Now  X  admits  a  singular  value  decomposition 

x  =  USV^ 

where  A  is  the  matrix  containing  the  singular  values  and  U,  V  are  orthogonal  matrices.  It 


A  AAA+AA+A 

c  =  USV  VS  V 

A  A  A  ^  A  + 

=  USS  V 


follows  that 


Next  consider  the  basis  function  vector  x  =  jx/t);  j  =  1,2,- ••.nj  which  can  be  written  as 

A  A  +  A 

X=U  X 

We  can  easily  show  that 

A  +  AA+A 

XX  =  U  XX  U 

A  A 

=  SS  =A 

after  some  straightforward  manipulations. 

Additionally,  we  have 

ft  A  ={^Zxf(t)dt 

A 

=  trc 

Comments 

The  above  analysis  permits  us  (at  least  in  principle)  to  enhance  SNR  by  orthogonal 
decomposition  methods.  Although  succinctly  stated,  the  numerical  computations  need  to  be 
carried  out  on  a  parallel  machine  to  test  the  efficacy  of  the  proposed  method. 


ABSTRACT 


The  statistics  of  a  harmonic  signal  (coherent  coinponenr)  mixed  with  a  random  back- 
giound  (incoheient  component)  oi  a  specified  spcxtral  prome  (pon’cr  spectrum)  is  still  a 
pioblem  of  inteiest.  The  puip(),s(-''  (.)i  the  pr(:'s('nt  j)aptn*  is  to  study  the  second-order  inten- 
sit}  statistics  of  such  a  signal/l)ackgr()und  sitnatieu  using  the  generalized  Karhunen-Loeve 
expansion^ 
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SECOND  ORDER  INTENSITY  STATISTICS 
OF  A  COHERENT  SIGNAL 
IN  THE  PRESENCE  OF  A  RANDOM  BACKGROUND 


1.  INTRODUCTION 


In  a  classic  paper  written  se\'eral  years  alter  the  actual  work  (done  during  World  War 
II),  Kac  and  Siegert  [1]  determined  the  exact  hrst-orcLer  statistics  (i.e..  probability  density 
and  moments)  of  a  scpiare  law  detector  for  a  Gaussian  landom  field,  and  foi  an  haimonic 
signal  buried  in  a  Gaussian  random  field.  Their  approach  involves  the  construction  of 
a  homogeneous  integral  equation  whose'  ke'riad  is  tlic'  coi'arianct'  lunction  of  the  landom 
field  using  what  is  now  termed  the  Karhuium-Loeve  expansion  [2.  3j.  although  we  can  also 
term  the  construction  the  Kac-Siegert  expansion  in  as  much  as  the}  deii\ed  it  indepen¬ 
dently.  The  eigenvalues  determiiK'  tlu'  prol)alhlity  density  function  of  the  detected  intensity 
(square  of  the  field  amplitude),  while  both  eigenvalues  and  eigenfunctions  are  needed  lOi 
the  probabilit}'  densit}'  of  the  int(‘usity  oi  the  sigmil  and  held.  Emerson  GJ  also  dibcusses 
the  problem  from  an  alternati\'('  viewpoint  using  th('  method  of  cumulants  to  avoid  solving 
the  associcited  homogeneous  int('gral  of  Kac  and  Siegert.  For  further  work  on  the  prob¬ 
lem.  see  Slepian  [5].  As  hlayer  and  kliddleton  [0]  ha^'c  pointed  out.  the  oiigmal  expansion 
method  of  Kac-Siegert  is  not  general  enough  to  handle  higher-order  statistics  such  as  prod¬ 
uct  moments.  Howe\'er.  Kac-Siegert  also  [)r('S('nr  another  method  of  solution,  the  direct 
method,  which  is  capable  of  handling  higher-order  statistics.  The  direct  method  lequires 
an  appropriate  transformation  to  ex})ress  the  output  in  terms  of  the  input;  the  statistics 
of  the  output  are  then  determined  b}'  suitable  additional  transformations  with  respect  to 
the  original  input  statistics,  Maym-  and  Middh'ton  exploit  this  approach  to  determine  the 
higher-order  statistics  of  the  output  due  to  a  square  law  detector  for  both  narrow-band 
and  broad-band  inputs.  Kac-Siegert  only  consider  the  broad-band  situation. 

The  purpose  of  the  present  paper  is  to  restudy  the  al^we  proldems  for  the  second-order 
intensity  statistics  using  a  generalisation  of  tlu‘  oi'iguial  Kac-Siegert  expansion  approach. 
I  emplc)V  two  point  detectors  to  iut('rrogat('  the  random  input  (with  and  without  the  hai- 
monic  signal  present).  T1k'S('  ck'tectors  op<'ratt‘  for  tlie  .banu'  tunc  int(*i\al  but  aict  dcla}cd 


relative  to  each  other  1\\’  a  varial)l(:'  time.  The  analysis  is  performed  \'ia  an  expansion  of  the 
random  field  over  these  two  disjoint  time  intervals  using  a  generalization  of  the  Karhunen- 
Loeve  series  developed  for  similar  problems  in  i)hoton  counting  [7-9]  and  lasei  speckle  [10]. 
In  this  approach,  a  homogeneous  integral  equation  is  constructed  over  the  two  disjoint  time 
intervals  wherein  the  two  detectors  operate.  Thc'  eigen\7i!ues  and  eigenfunctions  (which 
obey  an  unusual  orthogonality  condition)  are  e^•aluat(’d  and  used  to  fabricate  a  double 
generating  function:  from  this  doul)l('  geiK'rating  function  the  various  product  moments  of 
the  integrated  intensities  can  b("  obtaincal  In'  differentiation.  The  method  can  be  cairied 
out  for  the  third  order  and  higher  statistics.  se('  Blake  and  Barakat  ill]  for  the  thiid-oidei 
situation  in  the  context  of  photon  ('ountiiig.  Analysis  is  conrined  to  the  naiiow-band  situ¬ 
ation,  although  there  is  no  difficulty  in  studying  the  lu'oad-band  situation.  I  feel  that  the 
approach  via  the  generalized  I\.arhuiieii-Loev('  expansion  c.)fiers  a  nK)ie  ScUisfjing  physical 
picture  than  the  direct  nu.vhod  in  as  much  as  tlu'  detector  rime  intervals  and  delays  are 
an  inherent  part  of  the  analysis  \'ia  tlu'  associated  disjoint  integial  equation. 


2.  RANDOM  INPUT  FIELD 


The  coniplex- valued  amplitude  U[t)  of  the  total  disturbance  is  the  sum  of  a  determin 
istic  (or  coherent)  component  Uc(t)  and  a  random  l^ackground  term  Ub{t) 

U(t)  =  Uc[t) +  Ub{t)  .  (2.1) 


The  coherent  component  is  gi\-en  by 


Ucit)  =  <fce 


“  iu’e  t 


(2.2) 


where  <fc  const anr. 

The  random  background  UbiO  is  taken  ro  Ik*  a  /.ero-mean.  comjdex-valued,  spatially 
stationary  Gaussian  random  jn'oeess.  be.. 

Ubit]  —  \t)  +  iL'^  (t)  (2.3) 

where  is  the  stochastic  Hillrert  transiorui  ol  *(i).  Both  and  are  real¬ 
valued.  Now  and  have  the  same  Gaussian  probability  density  function 

(PDF)  and  are  statistical!}-  independent.  Furtliermore 

(h''’(t))  -(h‘>(t))  =0  (2.4) 

2 

(cFifab-’a.))  =  (Tli.iciPiO)  =  j  mti-t,)  (2.5) 

(c;"’ifiic;"u,))  =  o  (2.6) 
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where  a’l  is  the  variance  of  Ub{t)  and  rb[ti  -  h)  is  the  corresponding  correlation  function. 
0  <  |r6l  <  1.  Since  UbH)  is  a  Gaussian  random  process,  it  is  completely  characterized  by 
its  mean,  variance,  and  correlation  function. 

The  disturbance  U{t)  is  interrogated  by  two-point  detectors,  the  first  detector  op¬ 
erating  during  the  time  interval  {—T/2,T/2)  and  the  second  during  the  time  inteival 
(r  -  T/2,t  +  T/2)  where  t  is  the  variabb'  rime  delay.  The  integrated  intensities  Qj  are 

/.T/2 

O,  =  /  \U{tf  dt 

(2.7) 

/•r-tT/C 

()_,  =  /  \U(t)\-  dt  . 

Jt-T/1 

We  can  expand  Ubit)  in  a.  geiieiTdized  IvarhuiK'U-Loeve  series  ^6--9]  ovei  the  disjoint 
intervals 

=  (-r/2,r/2),  W  =  (T  -  T/2,r  +  T/2)  (2.S) 


SO  that 

CHO 

Ub[t)  =  Uki>k[t),  t  e  -4i  and  A.; 
=  0  ,  t  ^  Ai  and  A2 


(2.9) 


The  following  conditions  are  to  l)e  satisfied; 


1)  The  {Uk]  are  random  coefficients.  imU'iK'iident  of  t:  and  are  uncorrelated  Gaussian 
random  variables  (hence.  indep('iidenr  random  varial^les) 

{U,,U/)  =  albki  (2-10) 

where  {cr^}  are  as  vet  unknown,  real  noimegativc  constants.  It  is  essential  that  U{t)  be 
Gaussian,  for  if  it  is  not  then  tin-  exiiaiisioii  coefficients  {Uk}  vill  not  be  statistically 
independent,  although  they  will  still  be  uiicorrelated. 


J)  The  deterministic  functions  {'>pk(i)}  tt)  form  a  complete  orthonormal  set  ocei 


both  Ai  and  Ao- 


nonne^ 


The  precise  statement  of  the  (.)rtli()t!jouahty  condLitiou  is  cjuite  unusual.  Consider  the 
weighted  sum  of  the  integrated  inteiisitii's:  Ai  ill*’’  +  AjQo  '  '■''here  Ai.  A2  are  arbitrary  teal, 
leo'ative  parameters  appearing  in  the  two-fold  generating  function  of  the  integrated 


background  intensities  Q 


QblAijAj)  =  (expi -A|C‘i^'  -  AjO',  ’)) 


As  described  in  [i-O].  we  luu'c' 


(2.12) 


subject  to  the  requirement 


(j)  dt  =  81. 1 


(2.13) 


(b  =  X\ 


To  evaluate  the  unknown  constants  in  Eq.  (2.10).  we  must  construct  an  integial 
equation  whose  kernel  is  the  correlation  of  function  of  The  integial  equation  is 


'  .T--2  i>r+T/>\  /o-A-  ,  / 

Ai  J  +A2  J  jrbitx  -  hH'dM)  du  =  j  • 


(2.15) 


The  correlation  function  rb{t\  -  h)  <>f'  tin'  random  liackground  field  is  given  by 


■rt(f  1  -  i-2 )  = 


ti  -t-.-) 


9b{tl  -  t-2. 


(2.1G) 


G 


vhere  gb{ti  —  h)  is  the  eorrehitiou  fmicriou  ci'iireied  at  zero  frequency  and  oto  is  the 


frequency  at  which  the  power  spectrui 


n  lines 


diape)  of  the  background  radiation  is  a  max¬ 


imum.  If  we  set 


(2.17) 


then  Eq.  (2.15)  reads 


f»T  —  2  pr-TT/2\  /  ' 

Ai  /  -f-Aj  /  ]  gbiti  -  h)ov  U)  dt:  = 

J-T/l  Jt-T/>  I 


;2.1S) 


independent  of  a-u-  This  is  rhe  basic  inri'gral  ec,uati(m  for  dererinining  it  is  not  of 

the  standard  Fredholm  type  because  of  rlu'  presence  of  two  disjoint  domains  of  integration. 
A  second  difficulty  is  that  the  eigcnimlues  aj  are  implicit  iV.ncrions  of  Aj  and  A-.. 

The  two-fold  generating  fuiicrion  (/fc':Ai,Aj)  can  be  exprc'.ssed  as  an  infinite  product. 
To  prove  this  we  note  that  from  its  definition.  Etj.  i'2.11). 


Q6(  Ai ,  Aj  i 


f[{Uk}\c 


\  _  \.,  O'  ^  ' 


lld-Uk 


(2.19) 


where  d'Uk  =  and  /[{tT}]  the  joint  probability  density  function  of  the 

statistically  independent  Gaussian  r;nidoni  varialdes; 


/i{mi  =  ri  vc*' 


(2.20) 


Upon  substituting  Eq.  (2.15)  into  the  intc'grand.  we  encounter  a  standard  Gaussian  mte- 


o-ral.  It  can  Idc  shown  [6-9j  that 


QblAi,  A))  =  11 


2  1;  ) 


(2.21) 


Thus  far  we  have  only  considered  the  random  background  field  Ub{i).  In  the  next 
section  we  will  modify  the  analysis  to  include  tlu'  coherent  component  (signal). 


3.  RANDOM  INPUT  FIELD  WITH  COHERENT  SIGNAL 

The  starting  point  for  the  inclusion  of  the  signal  term,  Uc{t),  is  the  double  generating 
function  as  given  in  Eq.  (2.19)  vith  Q*'*  replaced  by  Q  ■.  In  particular 


AiQi  +  A2fl2  =  ^  \Ub(t)\''  dt  +  ^  'Tclt)!"  dt 


Ub[i)U:[t)dt-T 


u;{t)UAt)  dt . 


(3.1) 


(3.5) 


where 


6h.(Ai,A,) 


eft  . 


(3.6) 


Given  Eq.  (3.5).  we  have 

OO 

Q(Ai, A2)  =  exp{-|^c|'T'(Ai  +  A2)}  ]][  QfcfAnAo)  (3./) 

k^O 


where 


Qki^i }  ^2)  ~  i 


i__  n  (i-,a;-)ic4i-  ^-e,G:u,-uG,u:  ^2^^ 


iva 


k  j  -CO 


(3.S) 


To  evaluate  Qk-  we  note  tliat  since  Uk  =  14  +  il'14-  then 


CGlUk  --  icGkUt  =akV, 


(3.9) 


where 


ak  -  QGl  -r  icGk 

hk=^CcGl-icGk 


(3.10) 


Consequently 


Qr.= 


,-(14-a;'-)V^-a.n  /  ^-(l+cry  )Wi^-ibuW^  jy^k 


trtrr.  J- 


k  */-co 


(3.11) 


Both  integrals  are  known; 


/: 


G'‘ 


(1  +  ■) 


exp 


4(1 


(3.12) 


'‘CO  _o  *)  v^’tT 

g-(!-r-r)”A-  cos6a:B4  dWk  = 


\/il 


exp 


2  l2 

h  -  ^k 


[4(1  +  cr^ 


(3.13) 


0 


and  Eq.  (3.11)  reduces  to 


.  ,  1  \\^cGk\^4' 


(3.14) 


Thus  the  double  generating  function  for  the  background/signal  is 


Q(A„A,)  =  cxpi-liel-UA,  +  A,))  ■  Iltl  -  ^  11  (1+4)' 

k  =  ()  ^  =  0 

=  Q c  (  ^  i  J  ^'2  )  Qb (  ^  1  )  ) Qfcc  (  ^  I  1  ^'2  } 


where  Qc-  Qh  are  the  generating  functions  for  the  colierent.  Ijackground  components  re¬ 
spectively  and  Qic  the  generating  function  for  tlie  interaction  between  background  and 


The  product  moments  of  the  integrated  intensitie.s  can  be  obtained  by  differentiation 


of  <5(Ai ,  Ao): 


(3.16) 
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4.  DOUBLE  GENERATING  FUNCTION:  SMALL  APERTURE  REGIME 

It  should  be  obvious  that  the  larger  are  the  recei\-ing  apertures  the  greater  is  the 
smoothing  of  the  field  amplitudes,  Consc'tiuciitly  we  want  to  use  apertures  as  small  as  pos¬ 
sible,  and  now  assume  that  T  is  small  compared  to  th(:>  distance  over  which  the  correlation 
function  gb{t)  decays  to  its  e"'  value. 

Under  this  condition  only  the  k  =  0  terms  in  Qb  and  Qbc  contribute.  With  reference 
to  [7-9]  we  have 


rhA,,A2)  =  cr-TiA,  -rX,)-ra^T-\l  -  g-{r)]XiXj  . 


(4.1) 


This  leaves  only  the  Go  function  to  evaluate.  We  apply  the  mean  value  approximation  to 
evaluate  the  integrals  in  Eq.  (3.G).  therein'  obraimug 


Go  =  fc(A,r  + 


(4.2) 


where  (f>o{^)  =  <?o(")  =  constant  (=  b).  h'ote  that  |Gol'  is  a  quadratic  function  of  Ai  and 
A2.  Thus  we  have 


Q(Ai,A2)  =  (1  +  cr,])  '  cxp{-lTcl'T(Ai  -j- Asljexp 


|«^cGo|~<^o  ) 

(1  +  ^  J 


(4.3) 
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5.  PRODUCT  MOMENTS:  BACKGROUND/SIGNAL 

The  product  moments  of  the  integrated  intensities  were  obtained  by  differentiation 
of  <5(Ai,A9),  Eq.  (4.3),  according  to  Eq.  (3.16);  I  employed  the  symbol  naanipulation 
program,  SMP,  using  a  VAX  computer  to  eff'Ct  the  differentiations.  The  first  few  product 

moments  are: 

(Qj)  =  (Q.,)=cr--:sXl'  (5-1) 

(OjO.,)  =  -  g-{T)]  -r  •2ct-Tc1'  T  (5-2) 


(filO.d)  =  {Q,Ol)  =  lV-'II  +  2g-{r)] 

+  <7‘|,fcr'[4  +  ’2g-(T)  +  b  +  26co.s  Ar]  +  3cr’l,^cr‘  +  161*’  (5-3) 


=  4(7^  [1  +  45'(r )  +  g'l'’')) 

+  w  ^cos  At]  +  Sb] 

T  o-‘|6r‘[2  +  fif'l'r)  +  46  +  S6cos  At] 

+  4(7-161'’  + 161'  • 

Note  that  the  frequency  offset  A  only  appears  in  (QjQ^)  <tnd  {Q.1Q2) 

In  the  horaodyne  case  where  the  frcqnt'ncy  ol  the  signal  coincides  with  the  maximum 
of  the  power  spectrum  of  the  bt.ckground.  then  A  vanishes  and  the  product  moments  then 
depend  upon  gir)  only. 
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LDV  AND  TURBULENCE  ISSUES 


A  considerable  effort  during  this  contract  was  devoted  to  the  workings  of  LDV  itself  (as  a 
measurement  tool  for  turbulence).  The  basic  paper  in  this  area  is 

W.  George  and  J.  Lumley,  "The  laser-doppler  velocimeter  and  its  application  to  the 
measurement  of  turbulence",  J.  Fluid  Mech.,  60,  321-362  (1973). 

The  key  issue  is  the  measurement  of  the  doppler  shift  on  the  frequency  of  the  light  scattered  from 
small  particles  moving  in  the  turbulent  fluid.  George  and  Lumley  (as  well  as  others  since  their 
work)  have  concentrated  upon  this  doppler  shift  in  instantaneous  frequency  by  employing  ideas 
from  frequenc\-  modulation  (FM)  theory  due  to  Wang: 

J.  Lawson  and  G.  Uhlenbeck,  Threshold  Signals,  (McGraw-Hill.  New  York,  1950)  see 

chapter  13. 

w'ho  investigated  the  unmodulated  carrier  plus  noise  case.  Unfortunately  this  case  yields  a 
covariance  function  with  an  infinite  variance  (this  being  the  cause  of  all  the  troubles  in  the 
George-Lumley  approach). 

I  have  initiated  (and  essentially  carried  out)  the  analysis  when  the  Gaussian  noise  is 
broad-band  (as  in  the  turbulence  situation).  Unfortunately  the  analysis  is  much  more 
complicated  than  in  the  narrow-band  case  involving  many  multiple  integrals  over  Bessel 
functions  (which  reduce  to  simple  expressions  in  the  naiTOW-band  case).  The  mam  point, 
however,  is  that  in  this  case  the  covariance  function  has  a  finite  variance. 

Discussions  w  ith  Barbara  Sandler,  who  does  most  of  my  serious  programming,  have 
convinced  me  that  about  the  only  way  to  proceed  is  numerically.  We  have  developed  numerical 
integration  schemes  for  other  problems  that  can  be  used  for  the  broad-band  FM  situation.  Based 


on  these  considerations,  I  have  decided  to  postpone  writing  up  this  material  until  I  begin 
numerical  computations. 


